Upside Next Prefecture 分析

Section 1 データの読み込み・Function作成

・データフレームの説明 dt.weight : 漁獲量データ(1962~2015) dt.price : 産出額データ(2003~2015) translation_pref : 都道府県の翻訳データ tlanslation_fish : 魚種の翻訳データ multi_stock : 系群データ

#有効数字を999桁に設定
options(scipen=999)
#ディレクトリの変更
setwd("~/Google ドライブ/GD_Research_Upside/Upside_project/Data_Working")
#漁獲量データを読み込み、dt.weightに格納 
dt.weight<-read.csv("dataKaimenKenbetsu_modified_1962_2015.csv",stringsAsFactors = F)
#産出額データを読み込み、dt.priceに格納
dt.price<-read.csv("dataKaimenSanshutsugaku_modified.csv",stringsAsFactors = F)
#都道府県の翻訳データを読み込み、translation_prefに格納 
translation_pref<-read.csv("translation_prefecture.csv",stringsAsFactors = F)
#魚種の翻訳データを読み込み、translation_fishに格納  
translation_fish<-read.csv("translation_fish.csv",stringsAsFactors = F)
#系群データを読み込み、multi_stockに格納
multi_stock<-read.csv("multi_Stock_complete_All.csv", stringsAsFactors = F)
#佳奈恵さんのRunしたUpsideの結果データの読み込み
Temp_Data_Catch=read.csv("japan_20190118_mo_opt.csv", header=T,stringsAsFactors = F)
#徳永佳奈恵さんのデータで用いられている魚種名と行政データの魚種名統一用ファイルの読み込み
translation_kanae_kei<-read.csv("kanae_kei_fishlist.csv",stringsAsFactors = F)
#大海区は１つの県につき１つでは無いため、別ファイルで用意し、適宜使用
add_kaiku<-read.csv("translation_prefecture_region.csv",stringsAsFactors = F)#iranai

#小数点以下の数字の有効数字を調節するfunctionを作成(cutcut)
cutcut <- function(x,digits){
  return (floor(x * 10^digits) /10^digits)
}

Section 3 各系群の各県漁獲割合の算出

#期間は開始：2001年　終了　2010年　震災を含まないデータをつくるため。
Start_Year=2001
End_Year=2010

dt.comp <- dt.merge %>%
  subset(year>(Start_Year-1) & year<(End_Year+1)) %>%  ## 開始年と終了年の設定
  select(year,prefecture,alphabet_pref=alphabet,fish_stock,weight_t,value_million_JPY,Japanese,region_name) %>%
  group_by(fish_stock) %>%
  mutate(total_h = sum(weight_t)) %>%  ## 指定期間内の魚種別の総漁獲量
  mutate(total_p = sum(value_million_JPY,na.rm = T)) %>%  ## 指定期間内の魚種別の総産出額
  ungroup() %>%
  group_by(prefecture,fish_stock) %>%
  mutate(total_h_pref=sum(weight_t)) %>%  ## 指定期間内の県別かつ魚種別の総漁獲量
  mutate(total_p_pref=sum(value_million_JPY,na.rm = T)) %>%  ## 指定期間内の県別かつ魚種別の総産出額
  ungroup() %>%
  mutate(percentage_h = total_h_pref/total_h) %>%  ## 指定期間内の都道府県漁獲割合
  mutate(percentage_p = total_p_pref/total_p) %>%  ## 指定期間内の都道府県産出額割合
  unique()


DT::datatable(dt.comp,rownames = FALSE,filter = 'top')

Section 4 Upsideデータと過去のデータの統合

# 漁獲・生産額割合のみを抽出
dt.past <- dt.comp %>% 
  select(prefecture,alphabet_pref,fish_stock,percentage_h,percentage_p,region_name) %>%
  unique()

# 統合
dt.comp1 <- Temp_Data_Catch1 %>%
  left_join(dt.past,by="fish_stock") %>%
  mutate(pref_harvest = harvest1*percentage_h) %>% #pref_harvest = その年のその県の漁獲量：単位はt
  mutate(pref_profit = profit1*percentage_h) #pref_profit = その年のその県の生産額 ：単位は$

Section 5 Upside県別データ

今回からはHarvestとProfitを分けない。Profitを求める場合は以下の手順を踏むこと １dt.graphにtotal_pref_hが含まれている部分をtotal_pref_pに置き換える ２dt.graph2作成時のspread()でvalueをtotal_pref_pに置き換える ３時系列グラフ(下側)を作成する際にylim()を変更する

translation_pref_rm_na <- translation_pref %>% filter(pref_code>0)

dt.comp2 <- dt.comp1 %>%
  group_by(alphabet_pref,management,year) %>%
  mutate(total_pref_h = (sum(pref_harvest))/1000) %>% #その年、その県の総漁獲量(単位は1000t)
  mutate(total_pref_p = ((sum(pref_profit)/1000000))*110) %>% #その年の総産出額(単位は100万円：$1=¥110で計算)
  ungroup()

for (n in 1:nrow(translation_pref_rm_na)) {
  comp_name=translation_pref_rm_na[n,2]
  
  dt.graph <- dt.comp2 %>%
    filter(alphabet_pref==comp_name) %>%
    select(alphabet_pref,management,year,total_pref_h) %>% #産出額の場合はtotal_pref_p
    unique()
  
#comp_name1:2にグラフタイトルを格納
comp_name1 <- paste(comp_name,"Comparison of SQ and FMSY",sep="_")
comp_name2 <- paste(comp_name,"Timeseries",sep = "_")

#漁獲量においてFMSYがSQに追いつく年数の計算
dt.graph2 <- dt.graph %>%
  spread(key=management,value=total_pref_h) %>% #valueをtotal_pref_pにすると産出額になる
  mutate(Rate = ifelse(FMSY>SQ,"above","below")) ## 産出額の場合は下のylimの変更も必要

#head()を用いて並べる
dt.graph3 <- head(dt.graph2[dt.graph2$FMSY>dt.graph2$SQ,],1)
dt.graph4 <- if(nrow(dt.graph3)>0){if(dt.graph2$alphabet_pref=="hokkaido"){dt.graph3}else{rbind(dt.graph3,dt.graph4)}}else{if(dt.graph2$alphabet_pref=="hokkaido"){dt.graph2[1,]}else{rbind(dt.graph2[1,],dt.graph4)}}

# Write a bar graph econOpt / SQ
#dt.comp2_r$rateにeconOpt/SQ-1を格納。-1によってSQを基準(=1)として大小を正負で表すことができる
dt.graph2$rate <- (dt.graph2$FMSY/dt.graph2$SQ)-1

#グラフのy軸を最大値・最小値に合わせるために計算し、それぞれmax_r,min_rに格納
max_r <- max(dt.graph2$rate)
min_r <- min(dt.graph2$rate)

#値がマイナスの時には四捨五入すると軸が値よりも小さくなる場合があるため、マイナスの場合は切り上げを行う
if(min_r>0){min_r <- 0}else{min_r <- cutcut(min_r,digits=2)}

#fmsy/sqの棒グラフ
g1 <- ggplot(dt.graph2, aes(x = year,y = rate, fill = Rate)) +
  theme_light(base_family = "HiraKakuPro-W3")+ 
  labs(x="年", y="割合") +
  labs(title = comp_name1) +
  #scale_x_continuous(breaks = seq(2016,2065,7)) + #データが49個なので最大の素数7で分割する
  scale_x_continuous(breaks = c(2016,2026,2036,2046,2056,2065)) + #基本的に10刻みで最後だけ9個
  ylim(min_r,max_r) +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15,face="bold")) +
  scale_fill_manual(values=c("above"="#00552e","below"="#68be8d")) +
  geom_bar(stat = "identity",position = "dodge") +
  theme(legend.position = 'none')

print(g1)

#時系列データ(2016~2050)
# Extraction of maximum and minimum values (グラフ縦軸調整のための最大値と最小値の抽出)
max_h <- max(dt.graph$total_pref_h)
min_h <- min(dt.graph$total_pref_h)
min_h <- as.integer(floor(min_h)) #切り上げ
reach_year <- as.integer(dt.graph3[,2]) #FmsyがSQに追いつく年

data1<-ggplot(data = dt.graph,mapping=aes(x=year,y=total_pref_h,group=factor(management),colour=factor(management),shape = management)) +
  theme_light(base_family = "HiraKakuPro-W3") + 
  xlab("年") +
  ylab("漁獲量(1000t)") +
  labs(colour = "Management") +
  labs(title = comp_name2) +
  scale_x_continuous(breaks = seq(2016,2065,10)) +
  geom_vline(xintercept=reach_year,colour="red") +
  ylim(0,max_h) +
  #ylim(min_h,max_h) + #産出額の場合マイナスの値があるのでこちらを用いる
  geom_line() +
  scale_color_manual(values=c("FMSY"="#35a16b","SQ"="#ff2800","econOpt"="#0041ff")) +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18,face="bold")) +
  geom_point() +
  theme(legend.position = 'none')

print(data1)
}

Section 6 Upside地図

地図作成に用いるchoroplethrは拡張機能が乏しく、任意の範囲・区間での色分け等が行えないので、任意の範囲・区間のものに新たな値を付与し、その値ごとに色分けを行う必要がある。 (参照データフレーム ; dt.rate_h)

for (n in 1:nrow(translation_pref_rm_na)) {
  comp_name=translation_pref_rm_na[n,2]
  
  dt.map <- dt.comp2 %>%
    filter(alphabet_pref==comp_name) %>%
    filter(management %in% c("SQ","econOpt","FMSY")) %>%
    select(year,alphabet_pref,management,total_pref_h) %>% #profitの場合はtotal_pref_p
    unique()
  
  #最適な政策をOptim列に追加
  dt.map1 <- dt.map %>%
    spread(key=management,value=total_pref_h) %>% #profitの場合はtotal_pref_p
    mutate(Optim = ifelse(SQ>FMSY & SQ>econOpt,"SQ","Other")) %>%
    mutate(Optim = ifelse(FMSY>SQ & FMSY>econOpt,"FMSY",Optim)) %>%
    mutate(Optim = ifelse(econOpt>SQ & econOpt>FMSY,"econOpt",Optim))
  
  dt.map2 <- if(comp_name=="hokkaido"){dt.map1}else{rbind(dt.map2,dt.map1)}
}
  
#50年間合計での最適政策
dt.map3 <- dt.comp2 %>%
  select(year,alphabet_pref,management,total_pref_h) %>%
  unique() %>%
  group_by(alphabet_pref,management) %>%
  mutate(total_pref_h_mg = sum(total_pref_h)) %>% #managementごとの指定期間内総漁獲量
  ungroup() %>%
  select(alphabet_pref,management,total_pref_h_mg) %>%
  unique() %>%
  spread(key = management,value=total_pref_h_mg) %>%
  mutate(Optim = ifelse(SQ>FMSY & SQ>econOpt,"SQ","Other")) %>%
  mutate(Optim = ifelse(FMSY>SQ & FMSY>econOpt,"FMSY",Optim)) %>%
  mutate(Optim = ifelse(econOpt>SQ & econOpt>FMSY,"econOpt",Optim))

#choroplethrに対応するdata.frameを作成:50年間合計
dt.map4 <- select(.data=dt.map3,region=alphabet_pref, value = Optim)
#choropleth map (階級区分図) の作成
admin1_choropleth(country.name = "japan",
                  df           = dt.map4,
                  title        = "Optimal_harvest_total_50years",
                  legend       = "management",
                  num_colors   = ) +
　scale_fill_manual(values=c("FMSY"="#35a16b","SQ"="#ff2800","econOpt"="#0041ff"))

#choroplethrに対応するdata.frameを作成:50年後
dt.map4_1 <- dt.map2 %>%
  filter(year==2065) %>%
  select(region=alphabet_pref, value = Optim)
#choropleth map (階級区分図) の作成
admin1_choropleth(country.name = "japan",
                  df           = dt.map4_1,
                  title        = "Optimal_harvest_after_50years",
                  legend       = "management",
                  num_colors   = ) +
　scale_fill_manual(values=c("FMSY"="#35a16b","SQ"="#ff2800","econOpt"="#0041ff"))

Section 7 大海区ごとのデータ

0-県が複数の海区を跨ぐ　1-hokkaido_sea_of_japan_north_and_hokkaido_pacific_north 2-pacific_north 3-pacific_central 4-pacific_south 5-east_china_sea 6-sea_of_japan_west 7-sea_of_japan_north 8-seto_inland_sea

dt.kaiku <- dt.comp2 %>%
  left_join(add_kaiku,by=c("alphabet_pref","prefecture")) %>%
  select(year,alphabet_pref,fish_stock,management,region_kaiku_name,total_pref_h,total_pref_p,rate) %>%
  filter(management %in% c("SQ","econOpt","FMSY"))

#大海区ごとの総漁獲量・総産出額の計算。大海区を跨いでいる県は1/2して両方の海区に均等に配分
dt.kaiku1 <- dt.kaiku %>%
  mutate(total_pref_h = ifelse(rate==0.5,total_pref_h/2,total_pref_h)) %>% #pref_profit
  select(year,alphabet_pref,management,total_pref_h,region_kaiku_name) %>% #total_pref_p
  unique() %>%
  group_by(year,management,region_kaiku_name) %>%
  mutate(total_kaiku_h = sum(total_pref_h)) %>% #total_kaiku_p  total_pref_p
  ungroup() %>%
  select(year,management,region_kaiku_name,total_kaiku_h) %>%
  unique() 

kaiku_list <- dt.kaiku1 %>% select(region_kaiku_name) %>% unique()

for (n in 1:nrow(kaiku_list)) {
  comp_name=as.character(kaiku_list[n,1])

 dt.kaiku2 <- dt.kaiku1 %>%
   filter(region_kaiku_name==comp_name)
#時系列データ(2016~2050)
# Extraction of maximum and minimum values (グラフ縦軸調整のための最大値と最小値の抽出)
max_h <- max(dt.kaiku2$total_kaiku_h)
min_h <- min(dt.kaiku2$total_kaiku_h)
min_h <- as.integer(floor(min_h)) #切り上げ
comp_name1 <- paste("Timeseries_of_",comp_name,"_Harverst",sep="")

data1<-ggplot(data = dt.kaiku2,mapping=aes(x=year,y=total_kaiku_h,group=factor(management),colour=factor(management),shape = management)) +
  theme_light(base_family = "HiraKakuPro-W3") + 
  xlab("年") +
  ylab("漁獲量(1000t)") +
  labs(colour = "Management") +
  labs(title = comp_name1) +
  scale_x_continuous(breaks = seq(2016,2065,10)) +
  ylim(0,max_h) +
  #ylim(min_h,max_h) + #産出額の場合マイナスの値があるのでこちらを用いる
  geom_line() +
  scale_color_manual(values=c("FMSY"="#35a16b","SQ"="#ff2800","econOpt"="#0041ff")) +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18,face="bold")) +
  geom_point() +
  theme(legend.position = 'none')

print(data1)
}